home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / PV.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  12KB  |  500 lines
  1. /* 
  2.  * pv.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * P(oly_)V(ert)
  11.  */
  12.  
  13. #include <math.h>
  14. #include <malloc.h>
  15. #include <stdio.h>
  16. #include <stdlib.h>
  17. #include "ip.h"
  18.  
  19.  
  20. #define    ENLARGE        1              /* used for testing with xph.c */
  21. #define LOOP_COMPLETE    8.0
  22.  
  23. #define SQ2        1.414213562
  24. #define ZERO        0.0
  25.  
  26. #define ON        1
  27. #define OFF        0
  28.  
  29. #undef    DEBUG_CURV
  30. #undef    SHOW_PHIK
  31. #undef    SHOW_ARCL
  32.  
  33. #undef    ECHO_INPUT
  34. #undef    ECHO_VERT
  35. #define FABS(a)        ((a) > ZERO ? (a) : - (a))
  36.  
  37.  
  38.  
  39.  
  40. /*
  41.  * angular_bend()
  42.  *   DESCRIPTION:
  43.  *     angular_bend computes ZAHN's angular bend function
  44.  *     phi_star = phi_star(arc_length) for a given contour
  45.  *   ARGUMENTS:
  46.  *     d_phik(float *) change in direction as measured in PI/4
  47.  *     phi_k(float *) pointer to the differential values in angular direction
  48.  *       along the contour from an arbitrarily chosen starting point
  49.  *     nv(long) number of contour segments
  50.  *   RETURN VALUE:
  51.  *     total angular bend for the contour 
  52.  */
  53.  
  54. double
  55. angular_bend (float *d_phik, float *phi_k, long nv)
  56. {
  57.   int k, n_vertex = (int) nv + 1;
  58.   double tot_phi;
  59.  
  60. /*
  61.  * regarding turn angle at first vertex:
  62.  * first vertex is counted as i = 0 and i = nv;
  63.  */
  64.   *(phi_k + 0) = (float) 0.0;
  65.   for (k = 1; k <= nv; k++)
  66.     *(phi_k + k) = *(phi_k + k - 1) + *(d_phik + k - 1);
  67.  
  68.  
  69.   printf ("\n...total angle = %lf\n", tot_phi = *(phi_k + nv));
  70.  
  71.  
  72. #ifdef SHOW_PHIK
  73.   printf ("\n...phi_k:\n");
  74.   for (k = 0; k <= nv; k++)
  75.     printf ("%d %12.6f\n", k, *(phi_k + k));
  76. #endif
  77.  
  78.   if ((tot_phi != LOOP_COMPLETE) && (tot_phi != -LOOP_COMPLETE)) {
  79.     tot_phi = 99.0;
  80.     printf ("\n...total angular bend not equal to two PI! ");
  81.     printf ("-->something wrong!\n");
  82.   }
  83.   return (tot_phi);
  84. }
  85.  
  86.  
  87. /*
  88.  * arc_length()
  89.  *   DESCRIPTION:
  90.  *     arc_length computes normalized contour segment ("edge") lengths
  91.  *   ARGUMENTS:
  92.  *     d_lk(float *) pointer to contour length segments
  93.  *     arc_len(float *) pointer to normalized contour length segments to be filled
  94.  *     nv(long) number of contour segments
  95.  *   RETURN VALUE:
  96.  *     contour length of polygon
  97.  */
  98. double
  99. arc_length (float *d_lk, float *arc_len, long nv)
  100. {
  101.   int k, n_vertex = (int) nv + 1;
  102.   double c_len;
  103.  
  104. /*
  105.  * regarding turn length of edge ending in first vertex:
  106.  * first vertex is counted as i = 0 and i = nv;
  107.  */
  108.   *(arc_len + 0) = (float) 0.0;
  109.   for (k = 1; k <= nv; k++)
  110.     *(arc_len + k) = *(arc_len + k - 1) + (float) *(d_lk + k - 1);
  111.  
  112. /*
  113.  * normalize arc_length 
  114.  */
  115.   printf ("\n...contour_length = %lf\n", c_len = *(arc_len + nv));
  116.   for (k = 0; k <= nv; k++)
  117.     *(arc_len + k) *= (float) (100.0 / c_len);
  118.  
  119.  
  120. #ifdef SHOW_ARCL
  121.   printf ("\n...arc_len:\n");
  122.   for (k = 0; k <= nv; k++)
  123.     printf ("%d %12.6f\n", k, *(arc_len + k));
  124. #endif
  125.  
  126. /*
  127.  * before returning, correct contour_length for factor of two incurred
  128.  * in ZAHN's half integer representation of curvature points
  129.  * (see 'descript.c')
  130.  */
  131.   return (0.5 * c_len);
  132. }
  133.  
  134. /*
  135.  * vert_to_clen()
  136.  *   DESCRIPTION:
  137.  *     vert_to_clen computes arc_length (partial sums of segments lengths) for a polygon
  138.  *     defined in terms of the sequence of its vertices
  139.  *   ARGUMENTS:
  140.  *     v(point *) pointer to point struct (see ph.h)
  141.  *     s(float *) pointer to partial sums array (to be filled)
  142.  *     nv(long) number of vertices
  143.  *   RETURN VALUE:
  144.  *     contour length of polygon
  145.  */
  146.  
  147. double
  148. vert_to_clen (struct spoint *v, float *s, long nv)
  149. {
  150.   int i;
  151.   double delx, dely;
  152.   double c_len;
  153.  
  154. /*
  155.  * evaluate partial sums of edge lengths and accumulate in s
  156.  * (note: s de facto becomes a unit offset array)
  157.  */
  158.   *(s + 0) = (float) 0.0;
  159.   for (i = 0; i < nv; i++) {
  160.     delx = (double) ((v + i + 1)->x - (v + i)->x);
  161.     dely = (double) ((v + i + 1)->y - (v + i)->y);
  162.  
  163.     if (i == (nv - 1)) {        /* handle bdy closure */
  164.       delx = (double) ((v + 0)->x - (v + i)->x);
  165.       dely = (double) ((v + 0)->y - (v + i)->y);
  166.     }
  167.     if (FABS (delx) < 1.0)
  168.       delx = 0.0;
  169.     if (FABS (dely) < 1.0)
  170.       dely = 0.0;
  171.  
  172.     *(s + i + 1) = *(s + i) + (float) sqrt (delx * delx + dely * dely);
  173.   }
  174.  
  175. /* 
  176.  * normalize
  177.  */
  178.   c_len = *(s + nv);
  179.   for (i = 1; i <= nv; i++)
  180.     *(s + i) *= (float) (100.0 / c_len);
  181.  
  182.  
  183.  
  184. #ifdef DEBUG_LENGTH
  185.   printf ("\nVERT_TO_CLEN:\n");
  186.   printf ("...contour length of polygon: %lf\n", c_len);
  187.   printf ("...normalized cum. edge lengths:\n");
  188.   for (i = 0; i <= nv; i++)
  189.     printf ("s[%d]: %f\n", i, *(s + i));
  190. #endif
  191.  
  192.   return (c_len);
  193. }
  194.  
  195. /*
  196.  * curvature_energy()
  197.  *   DESCRIPTION:
  198.  *     curvature_energy computes curvature energy for contour
  199.  *     (ref: I.T.Young, J.E.Walker and J.E.Brown, 
  200.  *     Inform. and Control 25, 357-370(1974)
  201.  *   ARGUMENTS:
  202.  *     d_phik(float *) change in direction as measured in PI/4
  203.  *     d_lk(float *) pointer to contour length segments
  204.  *     nv(long) number of contour segments
  205.  *   RETURN VALUE:
  206.  *     contour length of polygon
  207.  */
  208.  
  209. double
  210. curvature_energy (float *d_phik, float *d_lk, long nv)
  211. {
  212.   register int k;
  213.   int del_phik, del_phikk;
  214.   float zahn_c_len;
  215.   float del_lk, del_lkk, del_pk;
  216.   float del_ce;
  217.  
  218.   double curv_en;
  219.  
  220. /*
  221.  *  evaluate curvature at ZAHN curvature points
  222.  */
  223.   zahn_c_len = (float) 0.0;
  224.   curv_en = 0.0;
  225.   del_lk = del_lkk = (float) 0.0;
  226.   del_phik = del_phikk = 0;
  227.   del_pk = (float) 0.0;
  228.  
  229. #ifdef DEBUG_CURV
  230.   printf ("\n...curvature_energy, extracted from ZAHN code:\n");
  231.   printf ("\n k    del_phik    del_lk     del_cesq    curv_en\n");
  232. #endif
  233.  
  234.   for (k = 1; k <= nv; k++) {
  235.     zahn_c_len += *(d_lk + k - 1);
  236.  
  237.     del_phik = (int) (*(d_phik + k - 1));
  238.     if (del_phik % 2 != 0)
  239.       del_lk = (float) (0.5 * SQ2);
  240.     else
  241.       del_lk = (float) 0.5;
  242.  
  243.     del_phikk = (int) (*(d_phik + k));
  244.     if (del_phikk % 2 != 0)
  245.       del_lkk = (float) (0.5 * SQ2);
  246.     else
  247.       del_lkk = (float) 0.5;
  248.  
  249.     del_pk = del_lk + del_lkk;
  250.     del_ce = ((float) del_phik) / del_pk;
  251.     curv_en += del_ce * del_ce;
  252.  
  253. #ifdef DEBUG_CURV
  254.     printf ("%d  %d   %f  %f  %f\n",
  255.             k, del_phik, del_lk, del_ce * del_ce, curv_en);
  256. #endif
  257.  
  258.   }
  259.  
  260. /*
  261.  * normalize
  262.  */
  263.   curv_en *= (45.0 * 45.0) / (zahn_c_len * ((float) nv));
  264.  
  265.   return (curv_en);
  266. }
  267.  
  268. /*
  269.  * find_interior_pt()
  270.  *   DESCRIPTION:
  271.  *     find_interior_pt finds an interior point of polygon,
  272.  *     given as a sequence of vertices
  273.  *   ARGUMENTS:
  274.  *     v(point *) pointer to vertices (see ph.h)
  275.  *     nv(int) number of vertices
  276.  *     xi, yi(int *) interior point found
  277.  *     e_out(int *) pointer to edge out array
  278.  *     tot_phi(double) total angular bend (currently not used)
  279.  *   RETURN VALUE:
  280.  *     1 if interior point found, 0 otherwise
  281.  */
  282.  
  283. int
  284. find_interior_pt (struct spoint *v, int nv, int *xi, int *yi, int *e_out, double tot_phi)
  285. {
  286.   int segm_found;
  287.   register int k;
  288.  
  289.   k = 0;
  290.   segm_found = 0;
  291.   while ((segm_found == 0) && (k <= nv)) {
  292.     if (tot_phi < 0.0) {
  293.       if ((*(e_out + k) > *(e_out + k + 1)) &&
  294.           (*(e_out + k + 1) > *(e_out + k + 2))) {
  295.         segm_found = 1;
  296.         *yi = ((v + k)->y + (v + k + 2)->y) / 2;
  297.         *xi = ((v + k)->x + (v + k + 2)->x) / 2;
  298.       }
  299.     }
  300.     else {
  301.       if ((*(e_out + k + 2) > *(e_out + k + 1)) &&
  302.           (*(e_out + k + 1) > *(e_out + k))) {
  303.         segm_found = 1;
  304.         *yi = ((v + k)->y + (v + k + 2)->y) / 2;
  305.         *xi = ((v + k)->x + (v + k + 2)->x) / 2;
  306.       }
  307.     }
  308.     k++;
  309.   }
  310.   return (segm_found);
  311.  
  312. }
  313.  
  314.  
  315. /*
  316.  * reconstruct polygon
  317.  */
  318. /*
  319.  * reconstruct_poly()
  320.  *   DESCRIPTION:
  321.  *     reconstruct_poly reconstructs a polygon,
  322.  *     given a tanential representation of the object
  323.  *   ARGUMENTS:
  324.  *     d_phik(double *) change in direction as measured in PI/4
  325.  *     d_lk(double *) pointer to contour length segments
  326.  *     nv(int) number of vertices
  327.  *     v(point *) pointer to vertices to be filles (see ph.h)
  328.  *     porg(point *) pointer to origin of polygon (see ph.h)
  329.  *     e_outk(int *) edge out
  330.  *     tot_phi(double) total angular bend (currently not used)
  331.  *     imgIO(Image *) pointer to Image struct (see tiffimage.h)
  332.  *     value(int) value to use for drawing
  333.  *   RETURN VALUE:
  334.  *     area in pixels of reconstructed polygon
  335.  */
  336.  
  337. unsigned int
  338. reconstruct_poly (float *d_phik, float *d_lk,
  339.                   long nv, struct spoint *v, struct spoint *porg,
  340.                   int e_outk, double tot_phi,
  341.                   Image * imgIO, int value)
  342. {
  343.   unsigned int p_area;
  344.  
  345.   int deltalk;
  346.   int *d_l, *e_out;
  347.   int *xv, *yv;
  348.   int nvc;
  349.  
  350.   int xi, yi;
  351.   int col_index;
  352.   int status;
  353.  
  354.   register int k;
  355.  
  356.  
  357.   if ((d_l = (int *) calloc (((int) nv + 1), sizeof (int))) == NULL)
  358.       exitmess ("\n...mem allocation for d_l failed\n", 1);
  359.   if ((e_out = (int *) calloc (((int) nv + 1), sizeof (int))) == NULL)
  360.       exitmess ("\n...mem allocation for edge_out failed\n", 1);
  361.  
  362.   nvc = nv + 1;
  363.  
  364. /*
  365.  * reconstruct polygon in form {xk, yk}, 
  366.  * pick arbitrary starting point  
  367.  */
  368.   *(d_l + 0) = 0;
  369.   *(e_out + 0) = e_outk;
  370.   (v + 0)->x = porg->x;
  371.   (v + 0)->y = porg->y;
  372.  
  373.  
  374. #ifdef ECHO_INPUT
  375.   printf ("\n...input data:\n");
  376.   printf ("...first vertex: (%3d,%3d)\n", (v + 0)->x, (v + 0)->y);
  377.   printf ("...edge_out at first vertex: %d\n", e_outk);
  378.   printf ("...total turn angle: %f\n", tot_phi);
  379.   for (k = 0; k <= nv; k++)
  380.     printf (" %d      %f      %f\n", k, *(d_phik + k), *(d_lk + k));
  381. #endif
  382.  
  383.   for (k = 1; k <= nv; k++) {
  384.     if (e_outk % 2 != 0)
  385.       deltalk = (int) ((*(d_lk + k - 1) / SQ2) * ENLARGE + 0.5);
  386.     else
  387.       deltalk = (int) ((*(d_lk + k - 1)) * ENLARGE + 0.5);
  388.  
  389.     switch (e_outk) {
  390.     case 0:
  391.       (v + k)->x = (v + k - 1)->x + deltalk;
  392.       (v + k)->y = (v + k - 1)->y;
  393.       break;
  394.     case 4:
  395.       (v + k)->x = (v + k - 1)->x - deltalk;
  396.       (v + k)->y = (v + k - 1)->y;
  397.       break;
  398.     case 2:
  399.       (v + k)->x = (v + k - 1)->x;
  400.       (v + k)->y = (v + k - 1)->y - deltalk;
  401.       break;
  402.     case 6:
  403.       (v + k)->x = (v + k - 1)->x;
  404.       (v + k)->y = (v + k - 1)->y + deltalk;
  405.       break;
  406.     case 1:
  407.       (v + k)->x = (v + k - 1)->x + deltalk;
  408.       (v + k)->y = (v + k - 1)->y - deltalk;
  409.       break;
  410.     case 5:
  411.       (v + k)->x = (v + k - 1)->x - deltalk;
  412.       (v + k)->y = (v + k - 1)->y + deltalk;
  413.       break;
  414.     case 3:
  415.       (v + k)->x = (v + k - 1)->x - deltalk;
  416.       (v + k)->y = (v + k - 1)->y - deltalk;
  417.       break;
  418.     case 7:
  419.       (v + k)->x = (v + k - 1)->x + deltalk;
  420.       (v + k)->y = (v + k - 1)->y + deltalk;
  421.       break;
  422.     }
  423.  
  424.     e_outk = e_outk + (int) (*(d_phik + k - 1));
  425.     if (e_outk < 0)
  426.       e_outk = e_outk + 8;
  427.     else if (e_outk >= 8)
  428.       e_outk = e_outk % 8;
  429.  
  430.     *(e_out + k) = e_outk;
  431.     *(d_l + k) = deltalk;
  432.   }
  433.   *(d_l + 0) = *(d_l + nv);
  434.   *(e_out + 0) = *(e_out + nv);
  435.  
  436.   for (k = 0; k <= nv; k++) {
  437.     (v + k)->x = (v + k)->x / 2;
  438.     (v + k)->y = (v + k)->y / 2;
  439.     *(d_l + k) = *(d_l + k) / 2;
  440.   }
  441.  
  442. #ifdef ECHO_VERT
  443.   printf ("\n    k     x_vert     y_vert     d_lk    edge_outk\n\n");
  444.   for (k = 0; k <= nv; k++)
  445.     printf (" %4d      %4d       %4d       %4d         %4d\n",
  446.             k, (v + k)->x, (v + k)->y, *(d_l + k), *(e_out + k));
  447. #endif
  448.  
  449.   if (((v + nv)->x != (v + 0)->x) || ((v + nv)->y != (v + 0)->y))
  450.     exitmess ("\n...polygon is not closed! something wrong!\n", 1);
  451.  
  452.  
  453. /*
  454.  * set up graphics and fill polygon, including its edges
  455.  */
  456.  
  457.   if ((xv = (int *) calloc ((size_t) (nv + 1), sizeof (int))) == NULL)
  458.       exitmess ("\n...mem allocation for xv failed\n", 1);
  459.   if ((yv = (int *) calloc ((size_t) (nv + 1), sizeof (int))) == NULL)
  460.       exitmess ("\n...mem allocation for yv failed\n", 1);
  461.  
  462.   for (k = 0; k <= nv; k++) {
  463.     *(xv + k) = (v + k)->x;
  464.     *(yv + k) = (v + k)->y;
  465.   }
  466. /*
  467.  * Draw and fill the polygon
  468.  */
  469.   p_area = draw_polyfill (xv, yv, nvc, imgIO, value);
  470.  
  471. /*
  472.  * move cursor to point inside the polygon 
  473.  */
  474.   if ((status = find_interior_pt (v, (int) nv, &xi, &yi, e_out, tot_phi)) == 0)
  475.     exitmess ("...no convex vertex found!\n", 1);
  476.  
  477. #ifdef FOO
  478. /*
  479.  * invoke flood function and determine pixel count (area) 
  480.  */
  481.   col_index = 0;
  482.   p_area = flood_region (xi, yi, col_index);
  483. #endif
  484.  
  485. /* 
  486.  * overlay border on displayed polygon 
  487.  */
  488.   col_index = GRAY;
  489.   draw_poly (xv, yv, nvc, imgIO, col_index);
  490. /*
  491.  * deallocate memory 
  492.  */
  493.   free (d_l);
  494.   free (e_out);
  495.   free (xv);
  496.   free (yv);
  497.  
  498.   return (p_area);
  499. }
  500.